\documentclass[paper=a4, fontsize=11pt]{scrartcl} % A4 paper and 11pt font size

\usepackage[T1]{fontenc} % Use 8-bit encoding that has 256 glyphs
\usepackage{fourier} % Use the Adobe Utopia font for the document - comment this line to return to the LaTeX default
\usepackage[english]{babel} % English language/hyphenation
\usepackage{amsmath,amsfonts,amsthm} % Math packages
\usepackage[UTF8]{ctex}
\usepackage{xcolor}
\usepackage{listings}
\usepackage{tikz}
\usepackage{lipsum} % Used for inserting dummy 'Lorem ipsum' text into the template
\usepackage{clrscode}
\usepackage{sectsty} % Allows customizing section commands
\usepackage[framed,numbered,autolinebreaks,useliterate]{mcode}
\usetikzlibrary{graphs}
\usetikzlibrary{shapes.arrows}
\allsectionsfont{\centering \normalfont\scshape} % Make all sections centered, the default font and small caps
%\allsectionsfont{\indent\normalfont\scshape} % Make all sections centered, the default font and small caps
\usepackage{indentfirst}
\setlength{\parindent}{2em}
\usepackage{fancyhdr} % Custom headers and footers
\pagestyle{fancyplain} % Makes all pages in the document conform to the custom headers and footers
\fancyhead{} % No page header - if you want one, create it in the same way as the footers below
\fancyfoot[L]{} % Empty left footer
\fancyfoot[C]{} % Empty center footer
\fancyfoot[R]{\thepage} % Page numbering for right footer
\renewcommand{\headrulewidth}{0pt} % Remove header underlines
\renewcommand{\footrulewidth}{0pt} % Remove footer underlines
\setlength{\headheight}{13.6pt} % Customize the height of the header

\numberwithin{equation}{section} % Number equations within sections (i.e. 1.1, 1.2, 2.1, 2.2 instead of 1, 2, 3, 4)
\numberwithin{}{section} % Number figures within sections (i.e. 1.1, 1.2, 2.1, 2.2 instead of 1, 2, 3, 4)
\numberwithin{table}{section} % Number tables within sections (i.e. 1.1, 1.2, 2.1, 2.2 instead of 1, 2, 3, 4)

\setlength\parindent{0pt} % Removes all indentation from paragraphs - comment this line for an assignment with lots of text

\usepackage[top=2cm, bottom=2cm, left=2cm, right=2cm]{geometry}
\usepackage[linesnumbered,boxed,lined,ruled,vlined]{algorithm2e}
\usepackage{algorithmicx}
\usepackage{algpseudocode}

\usepackage{graphicx}
\usepackage{float}
%----------------------------------------------------------------------------------------
%	TITLE SECTION
%----------------------------------------------------------------------------------------

\newcommand{\horrule}[1]{\rule{\linewidth}{#1}} % Create horizontal rule command with 1 argument of height

\title{
\normalfont \normalsize
\textsc{中国科学院大学}\ \textsc{计算机与控制学院} \\ [25pt] % Your university, school and/or department name(s)
\horrule{0.5pt} \\[0.4cm] % Thin top horizontal rule
\huge 算法设计与分析第五次作业 \\ % The assignment title
\horrule{2pt} \\[0.5cm] % Thick bottom horizontal rule
}

\author{黎吉国} % Your name

\date{\normalsize\today} % Today's date or a custom date

\begin{document}

\maketitle % Print the title
\newpage
\section{Load balance}
\noindent
You have some different computers and jobs. For each job, it can only be done on one of two specified computers.
The load of a computer is the number of jobs which have been done on the computer. Give the number of jobs and two computer ID for each job. You task is to minimize the max load.

(hint: binary search)\\
\textbf{解：}\\
设有M个任务，N台电脑，构造流网络如下：\\
\begin{enumerate}
  \item Vertex:$s,t,p[N],q[M]$,其中p[N]对应N台电脑，p[M]对应M个任务。
  \item Edge:$s\to p[k],q[m]\to t,,k=1,2,3\ldots,N,m=1,2,3,\ldots,M$.\\
  $p[k]\to q[m]\text{(if 任务m可以由电脑k完成)}$。
  \item $Capacity:c(s,p[k])=L,\text{(L表示最大负载)}$\\
  $c(p[k],q[m])=1$\\
  $c(q[m],t)=1$
  \item Obj:$min(L)$
\end{enumerate}
上面的流网络中，L是$[1,M]$的整数值，而且要满足$|f|=M$的要求（正好把任务做完），同时还要是最小的。即有$|f|_{L}=M && |f|_{L-1}<M$。
\subsection{pseudo-code}
\begin{algorithm}
  \caption{load-balance}
      \KwData{jobs[M],computer[N]}
      \KwResult{ min(max(load)) }
      \Begin{
        $G(E,V,C)=constructGraph(jobs[M],computer[N])$\\
        $L=M,ceiling=M,floor=0$;\\
        \While{1}
        {
          $f=maxFlow(G)$;\\
          \If{$floor==ceiling-1$}
          {
            $G(E,V,C)=reset(G(E,V,C),L)$;\\
            $f=maxFlow(G)$;\\
            \Break
          }
          \If{$f>=M$}
          {
            $ceiling=L$;\\
            $L=(floor+L)/2$;\\
          }
          \Else
          {
            $floor=L$;\\
            $L=(L+ceiling)/2$;\\
          }
          $G(E,V,C)=reset(G(E,V,C),L)$;\\
        }
        \Return $f$;
      }
\end{algorithm}
\subsection{Correction}
在这个流网络中，汇点的特点使得最大流不会超过M，同时，每个电脑的负载不超过L，当最大流为M时，所有任务都被做完，当最大流小于M时，说明L对负载的限制使得一部分任务没有做完，
所以最小的保证最大流为M的L就是我们要找的最佳分配。
\subsection{Complexity}
计算最大流的复杂度为$O(V^2 E)=O((M+N)^2(M+N+2M))=O((M+N)^3)$，循环次数为$\log_2 M$，所以整体复杂度为$O(\lg M(M+N)^3)$.
\section{Matrix}
\noindent
For a matrix filled with $0$ and $1$, you know the sum of every row and column. You are asked to give such a matrix which satisfys the conditions.

\textbf{解：}\\
设矩阵的行和，列和分别为$Sum_r[M],Sum_c[N]$，满足$\sum_{k=1}^M Sum_r[k]=\sum_{m=1}^N Sum_r[m]$,可以构造如下流网络：
\begin{enumerate}
  \item Vertex:$s,t,V_r[M],V_c[N]$,其中$V_r[M]$对应M行，$V_c[N]$对应N列。
  \item Edge:$s\to V_r[k],V_c[m]\to t,V_r[k]\to V_c[m],k=1,2,3\ldots,M,m=1,2,3,\ldots,N$.
  \item Capacity:\\
  $c(s,V_r[k])=Sum_r[k]$\\
  $c(V_r[k],V_c[m])=1$\\
  $c(V_c[m],t)=1$
  \item 该网络的一个最大流如果等于$\sum_{k=1}^M Sum_r[k]$,则该最大流中$V_r[k]\to V_c[m]$的边便是满足要求的一个矩阵，否则，没有满足要求的矩阵。
\end{enumerate}
\subsection{pseudo-code}
\begin{algorithm}
  \caption{Matrix}
      \KwData{$Sum_r[M],Sum_c[N]$}
      \KwResult{ Matrix $A$ }
      \Begin{
        $G(E,V,C)=constructGraph(Sum_r[M],Sum_c[N])$\\
        $f=maxFlow(G)$;\\
        \If{$|f|=\sum_{k=1}^M Sum_r[k]$}
        {
          $A(k,m)=f(V_r[k],V_c[m]),k=1,2,3,\ldots,M,m=1,2,3,\dots,N$;\\
          \Return $A$;
        }
        \Else
        {
          \Return "Failed";
        }
      }
\end{algorithm}
\subsection{Correction}
当该网络的最大流为$\sum_{k=1}^M Sum_r[k]$时，有
\[
\begin{split}
  f(s,V_r[k])&=Sum_r[k], k=1,2,3,\ldots,M\\
  f(V_c[m],t)&=Sum_c[m], m=1,2,3,\ldots,N\\
\sum_{m=1}^N f(V_r[k],V_c[m])&= f(s,V_r[k])=Sum_r[k]\\
\sum_{k=1}^M f(V_c[k],V_c[m])&= f(V_c[m],t)=Sum_c[m]
\end{split}
\]
所以由$f(V_c[k],V_c[m])$组成的矩阵是一个满足要求的矩阵。
\subsection{Complexity}
给算法的复杂度和计算最大流的复杂度相等,为$O(V^2E)=O((M+N)^2MN)$

\section{unique cut}
\noindent
Let $G = (V, E)$ be a directed graph, with source $s \in V$ , sink $t \in V$ , and nonnegative edge capacities $c_e$. Give a polynomial-time algorithm to decide whether $G$ has a unique minimum $s − t$ cut.

\textbf{解：}\\
该问题，先考虑如何由一个最大流得到对应的最小割。得到最大流时，残余网络中没有从$s$到$t$的通路，也即是说，一部分节点可以由$s$出发到达，令一部分节点可以从它出发达到$t$，这两个集合是不相交的（一旦相交则会形成一条$s\to t $的通路）。
在我们做图的割时，只要保证这两个集合分别再两个子集中即可。但是如果存在节点，既不能由$s$到达，又不能有它出发达到$t$，那么在做割时就可以任意处理，这就造成了最大流所对应的最小割的不唯一性。
所以我们可以通过检测最大流对应的残余网络的节点的特性来判断最小割的唯一性。
\subsection{pseudo-code}
\begin{algorithm}
  \caption{Matrix}
      \KwData{$G(V,E,C)$}
      \KwResult{ "yes" or "no" }
      \Begin{
        $G_f(E,V,C)=maxFlow(G(E,V,C))$\\
        $S=\{ v|v\in V-\{ s,t \},\text{v is accessible from s in }G_f \}$;\\
        $T=\{ v|v\in V-\{ s,t \},\text{t is accessible from v in }G_f \}$;\\
        \If{$S+T == V$}
        {
          \Return "yes";
        }
        \Else
        {
          \Return "no";
        }
      }
\end{algorithm}
\subsection{Correction}
分析将上文。
\subsection{Complexity}
计算最大流的复杂度为$O(V^2E)$,计算两个子集的复杂度和最短距离的复杂度相同，是$O(VE)$，
所以整个算法的复杂度是$O(V^2E)$.

\section{Maximum flow}
Another way to formulate the maximum-flow problem as a linear program is via flow decomposition. Suppose we consider all (exponentially many) s-t paths $p$ in the network G, and let $f_p$ be the amount of flow on path $p$. Then maximum flow says to find

\[
\begin{array}{rrrrrr}
 \max & z & = &\sum{f_p}   & & \\
 s.t. &  \sum_{e \in p}f_p  & \leq & u_e ,& for \quad all \quad edge \quad e &   \\
      & f_p  & \geq & 0  & &  \\
\end{array} \nonumber
\]

(The first constraint says that the total flow on all paths through e must be less than ue.) Take the dual of this linear program and give an English explanation of the objective and constraints.

\textbf{解：}\\
\[
\begin{array}{rrrr}
  \min & \sum u_e y_e &\\
  s.t & \sum_{e\in P}y_e&\ge 1&\text{for all path }P\\
  &y_e&\ge 0&
\end{array}
\]
该线性规划求解最小割问题，$y$对应于图中的边，$y=1$表示该边是最小割的关建边，即该边连接两个子集。在每一条$s\to t$的通路中，至少有一个边被割，所以通路上$y$的和至少为1.该优化问题最终得到一个最小割。

\newpage
\section{Ford-Fulkerson algorithm}
Implement Ford-Fulkerson algorithm to find the maximum flow of the following network, and list your intermediate steps.
Use you implementation to solve problem 1 and show your answers.
\noindent
INPUT: ($N$, $M$) means number of jobs and computers. Next $N$ line, each line has two computer ID for a job. see more detial in the file problem1.data.

\noindent
OUTPUT: the minimum number of the max load.
\\
代码以及测试结果如下：
\lstset{language=Matlab}%代码语言使用的是matlab
\lstset{breaklines}%自动将长的代码行换行排版
\lstset{extendedchars=false}%解决代码跨页时，章节标题，页眉等汉字不显示的问题
\begin{lstlisting}[frame=single]
function rv=load_balance(filename);
rv=0;
%read data
fid=fopen(filename,'r');

while(1)
    %read txt
    while(1)
        if feof(fid)
            return;
        end
        line=fgetl(fid);
        if length(line)>1 && line(1)~='#'
            break;
        end
    end
    data_temp=sscanf(line,'%d %d');

    M=data_temp(1);%job
    N=data_temp(2);%computer
    job_com=zeros(M,2);%specified computer
    for idx=1:M
         line=fgetl(fid);
         if length(line)>1 && line(1)~='#'
            data_temp=sscanf(line,'%d %d');
            job_com(idx,1)=data_temp(1);
            job_com(idx,2)=data_temp(2);
         end
    end
    %construct a graph
    V_num=1+M+N+1;
    Capacity=zeros(V_num,V_num);%0:no edge

    Capacity(V_num-M:V_num-1,V_num)=1;
    for idx=1:M
        Capacity(1+job_com(idx,1),1+N+idx)=1;
        Capacity(1+job_com(idx,2),1+N+idx)=1;
    end
    %init
    L=1;
    data_floor=1;
    data_ceiling=M;
    while(1)
        Capacity(1,2:N+1)=L;
        [C_f f]=Ford_Fullerson(Capacity);

        %devide
        if sum(f(1,:),2)-sum(f(:,1),1)>=M
            data_ceiling=L;
            L=floor((L+data_floor)/2);
        else
            data_floor=L;
            L=floor((L+data_ceiling)/2);
        end
        %exist
         if data_ceiling-data_floor<1.5
             %disp(f);
             disp('min max load:');
             disp(data_ceiling);
             break;
         end
    end
end
end
\end{lstlisting}
其中最大流算法的实现如下：
\lstset{language=Matlab}%代码语言使用的是matlab
\lstset{breaklines}%自动将长的代码行换行排版
\lstset{extendedchars=false}%解决代码跨页时，章节标题，页眉等汉字不显示的问题
\begin{lstlisting}[frame=single]
  function [C_f f]=Ford_Fullerson(C);
   %[C_f f]=Ford_Fullerson(C)
   %C:capacity of the network
   %C_f:residual network
   %f:flow

  C_size=size(C);
  f=zeros(C_size);
  C_f=C;
  %
  while(1)
   %get a path
   P=find_path(C_f);%P=[u,v;]
  if isempty(P)
      break;
  end
  c_f=min(C_f(P(:,1)+(P(:,2)*C_size(2)-C_size(2))));%flow of P
  f_c=zeros(C_size);
  f_c(P(:,1)+(P(:,2)*C_size(2)-C_size(2)))=c_f;%cloumn
   %f=f+f'
   f=f+f_c;
   %residual graph
   C_f(P(:,1)+(P(:,2)*C_size(2)-C_size(2)))=C_f(P(:,1)+(P(:,2)*C_size(2)-C_size(2)))-c_f;
   C_f(P(:,2)+(P(:,1)*C_size(2)-C_size(2)))=C_f(P(:,2)+(P(:,1)*C_size(2)-C_size(2)))+c_f;
  end
  end

  function P=find_path(C);
  C_size=size(C);
  if C_size(1)<2 || C_size(1)~=C_size(2)
      P=0;
      return;
  end
  L=C_size(1);
  S=1;
  T=2:L;
  P_record=[];
  while 1
      [max_value1 max_idx1]=max(C(S,T),[],1);
      [max_value2 max_idx2]=max(max_value1,[],2);
      max_value=max_value2;
      max_idx=[max_idx1(max_idx2) max_idx2];
      max_path=[S(max_idx(1)) T(max_idx(2))];
      if max_value==0
          P=[];
          return;
      end
      %[min_value min_idx]=min(C.* (C(S,T)>0));
      P_record=[P_record;max_path];
      if max_path(2)==L
          break;
      end
      S=[S max_path(2)];
      T=[T(1:max_idx(2)-1),T(max_idx(2)+1:length(T))];
  end
  %back for path
  P=[];
  target=L;
  for idx=size(P_record,1):-1:1
      if P_record(idx,2)==target
          P=[P_record(idx,:); P];
          target=P_record(idx,1);
      end
      if target==1
          break;
      end
  end
  end
\end{lstlisting}
测试结果如下：

\lstset{language=Matlab}%代码语言使用的是matlab
\lstset{breaklines}%自动将长的代码行换行排版
\lstset{extendedchars=false}%解决代码跨页时，章节标题，页眉等汉字不显示的问题
\begin{lstlisting}[frame=single]
>> load_balance('problem1.data');

min max load:
     2

min max load:
     1

min max load:
     2

min max load:
     3

min max load:
     1

min max load:
     1

min max load:
     3

min max load:
     1

min max load:
     1

min max load:
     1
\end{lstlisting}

\section{Push-relabel}
\noindent
Implement push-relabel algorithm to find the maximum flow of a network, and list your intermediate steps.
Use your implementation to solve problem 2 and write a check problem to see if your answer is right.

\noindent
INPUT: Numbers of rows and columns. And the sum of them. See more detial in the file problem2.data.

\noindent
OUTPUT: The matrix you get. Any one satisfy the conditions will be accept.

\\
读取数据并使用Push\_Relabel算法求解最大流，实现如下：
\begin{lstlisting}[frame=single]
  function A=Matrix(filename);
  A=0;
  fid=fopen(filename,'r');
  while 1
      while 1
          if feof(fid)
              return;
          end
          line=fgetl(fid);
          if isempty(line)|| line(1)=='#'
              continue;
          else
              break;
          end
      end
      data_temp=sscanf(line,'%d%d');
      row =data_temp(1);
      col =data_temp(2);
      format=repmat('%d',1,row);
      line=fgetl(fid);
      Sum_row=sscanf(line,format);
      format=repmat('%d',1,col);
      line=fgetl(fid);
      Sum_col=sscanf(line,format);
      %construct a graph
      V=1+row+col+1;
      Capacity=zeros(V,V);
      Capacity(1,2:1+row)=Sum_row';
      Capacity(1+row+1:1+row+col,V)=Sum_col;
      Capacity(2:1+row,1+row+1:1+row+col)=1;
      [C_f f]=Push_Relabel(Capacity);
      A=f(2:1+row,1+row+1:1+row+col);
      disp(A);
  end
  end
\end{lstlisting}
其中，Push\_Relabel算法实现如下：
\begin{lstlisting}
function [C_f f]=Push_Relabel(C);
%[C_f f]=Push_Relabel(C);
%C:capacity
%C_f:residual network
%f:maxflow

C_size=size(C);
if C_size(1)==0 || C_size(1)~=C_size(2)
    C_f=[];
    f=[];
    return;
end
V=C_size(1);
V_data=zeros(V,2);
%index define
H=1;%height
E=2;%esceed
%init V_data
C_f=C;
V_data(1,H)=V;
V_data(1,E)=-sum(C_f(1,:),2);
V_data(2:end,E)=C_f(1,2:end)';
C_f(1,2:end)=C_f(1,2:end)-V_data(2:end,E)';
C_f(2:end,1)=C_f(2:end,1)+V_data(2:end,E);
while 1
    %find a vertex which is overflow
    V_idx=find(V_data(1:end-1,E)>0,1);
    if isempty(V_idx)
        break;
    end
    %push or relabel
    V_accessible=find(C_f(V_idx,:)>0);
    [min_H min_idx]=min(V_data(V_accessible,H));
    if min_H<V_data(V_idx,H)%push
        push_idx=V_accessible(min_idx);
        push_value=min(V_data(V_idx,E),C_f(V_idx,push_idx));
        V_data(push_idx,E)=V_data(push_idx,E)+push_value;
        V_data(V_idx,E)=V_data(V_idx,E)-push_value;
        C_f(V_idx,push_idx)=C_f(V_idx,push_idx)-push_value;
        C_f(push_idx,V_idx)=C_f(push_idx,V_idx)+push_value;
    else%relabel
        V_data(V_idx,H)=min_H+1;
    end
end
f=C-C_f;
f(f<0)=0;
end
\end{lstlisting}
测试过程及部分结果如下：
\begin{lstlisting}
  >> Matrix('problem2.data');

     1     0     0     0     0     0     1     1     1     1
     1     0     0     0     0     0     1     1     1     1
     1     0     0     0     1     1     1     1     1     1
     1     0     1     0     0     1     1     1     1     1
     0     1     1     1     1     1     0     0     0     1
     0     1     1     0     0     1     0     0     0     0
     0     1     1     1     1     0     0     0     0     1
     0     1     1     1     1     1     1     0     0     1
     1     1     1     1     1     1     1     0     0     0
     1     1     1     0     0     0     0     0     0     0
\end{lstlisting}

\end{document}
